Analyses of plasma metabolites using a high performance four-channel CIL LC-MS method and identification of metabolites associated with enteric methane emissions in beef cattle

Reducing enteric methane (one greenhouse gas) emissions from beef cattle not only can be beneficial in reducing global warming, but also improve efficiency of nutrient utilization in the production system. However, direct measurement of enteric methane emissions on individual cattle is difficult and expensive. The objective of this study was to detect plasma metabolites that are associated with enteric methane emissions in beef cattle. Average enteric methane emissions (CH4) per day (AVG_DAILYCH4) for each individual cattle were measured using the GreenFeed emission monitoring (GEM) unit system, and beef cattle with divergent AVG_DAILYCH4 from Angus (n = 10 for the low CH4 group and 9 for the high CH4 group), Charolais (n = 10 for low and 10 for = high), and Kinsella Composite (n = 10 for low and 10 for high) populations were used for plasma metabolite quantification and metabolite-CH4 association analyses. Blood samples of these cattle were collected near the end of the GEM system tests and a high performance four-channel chemical isotope labeling (CIL) liquid chromatography (LC) mass spectrometer (MS) method was applied to identify and quantify concentrations of metabolites. The four-channel CIL LC-MS method detected 4235 metabolites, of which 1105 were found to be significantly associated with AVG_DAILYCH4 by a t-test, while 1305 were significantly associated with AVG_DAILYCH4 by a regression analysis at p<0.05. Both the results of the t-test and regression analysis revealed that metabolites that were associated with enteric methane emissions in beef cattle were largely breed-specific whereas 4.29% to 6.39% CH4 associated metabolites were common across the three breed populations and 11.07% to 19.08% were common between two breed populations. Pathway analyses of the CH4 associated metabolites identified top enriched molecular processes for each breed population, including arginine and proline metabolism, arginine biosynthesis, butanoate metabolism, and glutathione metabolism for Angus; beta-alanine metabolism, pyruvate metabolism, glycolysis / gluconeogenesis, and citrate cycle (TCA cycle) for Charolais; phenylalanine, tyrosine and tryptophan biosynthesis, phenylalanine metabolism, arginine biosynthesis, and arginine and proline metabolism for Kinsella Composite. The detected CH4 associated metabolites and enriched molecular processes will help understand biological mechanisms of enteric methane emissions in beef cattle. The detected CH4 associated plasma metabolites will also provide valuable resources to further characterize the metabolites and verify their utility as biomarkers for selection of cattle with reduced methane emissions.


Introduction
Climate change and global warming have been a concern for all humans since 1970s, and the main cause, greenhouse gas emissions, has attracted enormous attention with discussions in political, environmental, technological, and cultural areas [1].Greenhouse gas mainly includes carbon dioxide (CO 2 ), enteric methane (CH 4 ), nitrous oxide (N 2 O), etc., in which the impact of CH 4 on the climate change is more than 25 times greater than CO 2 .According to a widely quoted figure from United Nations Food and Agriculture Organization (FAO), livestock, particularly cattle, are responsible for 14.5 percent of global human-induced greenhouse gas emissions [2].And due to an expected doubling demand of global milk and meat by 2050, CH 4 emissions from livestock is predicted to substantially increase [3].In Canada, it is estimated that beef cattle produce up to 20 million tons carbon dioxide equivalent (CO 2 eq) of methane per year as a result of enteric fermentation, accounting for the great influence on the greenhouse effect [4].
Cattle release 6% of their ingested energy as eructated methane [5], thus reducing methane emissions not only can be beneficial in reducing carbon foot print of the industry, but also improve nutrition utilizing efficiency and decrease production costs.The mechanism of methane emissions from cattle remains to be fully revealed, but it has been shown that cattle methane emission is in part controlled by animal genetics, offering an opportunity to reduce methane emissions through genetic selection.Conventional genetic selection of low methane emission beef cattle based on directly measured phenotypes is challenging and expensive because methane emission data collection is time-and labor-intensive [6].Marker assisted selection (MAS) or genomic selection is an alternative approach with identified predictive markers from DNA, transcriptomics, metabolomics or other types [7].
Of those, metabolomics is an emerging field and a powerful tool to characterize complex biochemical phenotypes.It has the potential to reveal promising biomarker candidates associated with biomolecular processes of methane emissions.The major issues of current metabolomics studies include low metabolome coverage as well as not high quantification accuracy.As a result, new techniques are required.Previously, we reported the use of high-performance chemical isotope labeling (CIL) LC-MS method for the profiling of amine/phenol submetabolome [8].This derivatization method can significantly improve LC separation, MS sensitivity, and also provide precise relative quantification result.The power and promising future of this method have been realized from its successful application in various samples and areas [9][10][11][12][13].In order to further improve metabolome coverage, recently, we have developed other CIL approaches targeting different submetabolomes, such as carboxylic acids [14], hydroxyls [15], ketones and aldehydes [16].By using this divide-and-conquer strategy and combining all these four-channel methods together, we could gain a coverage as high as 90% of potential whole metabolome [17].Thus, in this work, we applied the high performance four-channel CIL LC-MS method to detect and quantify plasma metabolites from cattle that were measured for enteric methane emissions.Subsequently, we conducted association analyses to identify metabolites and enriched molecular processes that are significantly associated with enteric methane emissions in beef cattle.

Animal populations and management
Animals used in this study were from three beef cattle populations including Angus, Charolais, and Kinsella Composite (KC) that are located at the Roy Berg Kinsella Ranch, University of Alberta.The three beef cattle populations were previously described [18,19].Briefly, the Angus and Charolais cow herds were maintained via breeding using registered bulls with their pedigree information maintained by the Canadian Angus or Charolais association, respectively.The KC cow herd were produced from crosses among 3 composite cattle lines, namely beef synthetic 1, beef synthetic 2, and dairy × beef synthetic (DBS).The beef synthetic 1 was composed of 33% Angus, 33% Charolais, and about 20% Galloway, with the remainder from other beef breeds.The beef synthetic 2 composite was made up of about 60% Hereford and 40% other beef breeds.The dairy × beef synthetic was composed of approximately 60% dairy breeds (Holstein, Brown Swiss, or Simmental) and 40% beef breeds, mainly Angus and Charolais.The KC cow herd was maintained via breeding using Angus, Charolais, or University of Alberta Hybrid bulls that were produced from the three beef synthetic lines.All the three cow herds were bred between June to August and their calves were born between March to May.Calves remained with their dams until they were weaned at 6-7 months of age.All animals were raised and managed following the Canadian Council of Animal Care (CCAC) guidelines on the care and use of farm animals in research teaching and testing [20], and all the experimental procedures applied to the animals were approved by the University of Alberta Animal Care and Use Committee (AUP00000777).

Methane emissions data collection
Methane emissions were measured on Angus and Charolais heifers after weaning and on KC mature cows using the GEM (GreenFeed emissions monitoring) system at the Roy Berg Kinsella Ranch along with a feedlot test trial.All the animals in the three breed population were fed ad libitum once daily (8am) with the mixed ration (as fed basis) during the feedlot test and during adaptation to diet (25 d before test) period.Ingredient composition of the mixed ration offered to animals in the Angus and Charolais populations was 71.8% silage, 19% oat grain, 4.5% corn dried distillers grains with solubles (cDDGS), and 4.7% Feedlot 30-Rumensin (F30).For the KC population, the ingredient composition was 85% Barley Silage, 10% Barley, 5% F30.The F30 contained crude Protein (min) 30.0%;Equivalent Crude Protein from Non-protein sources (max) 10.0%; Crude Fat (min) 2.5%; Crude Fibre (max) 8.0%; Flourine (max) 120 mg/kg; Sodium (act) 2,0%; Calcium (act) 8.5%; Phosphorus (act) 0.7%; Magnesium (act) 0.4%; Potassium (act) 0.9%; Sulfur (act) 0.4%; Iron (act) 135 mg/kg; iodine (act) 30 mg/kg; Copper (act) 310 mg/kg; Manganese (act) 840 mg/kg; Zine (act) 1,440 mg/kg; Cobalt (act) 6.0 mg/kg: Vitamin A (min) 100,000 LU./kg: Vitamin D (min) 10,000 I.U./kg:Vitamin E (min) 300 I.U./ kg and monensin 22mg/kg.Details of measuring methane emissions using the GEM system were described by Manafiazar et al. [21].Briefly, the GEM system dispenses feed pellets to attract the animals.When an individual animal visits GEM, the system reads the animal's radio-frequency identification tag, and the beginning, and end time of each visit are recorded.Any measurements with a minimum of 3 min per each visit time were considered as acceptable observations for monitoring emissions.Animals had free access to the GEM system 24 h/day while on drylot, although pellet intake was restricted such that each animal received a maximum of six drops per GEM visit, and pellets were dropped at 36 second intervals.After the sixth drop, animals were required to wait 4 h until their next six drops, resulting in a maximum of 36 drops per animal per day with six possible visits, and drop sizes averaged 35.4 g drop −1 (SD = 0.30).Pellets used in the GEM unit consisted of barley, beef vitamin-trace mineral premix, calcium carbonate, corn distiller's grains screenings, sodium chloride, wheat/wheat middlings, and zinc chelate (MasterFeeds Inc., Red Deer, AB, Canada).Continuous and negative air flow from a system fan draw air past the animal's nose and mouth when it enters the shroud, thus mixing air with respired and (or) eructated CH 4 and CO 2 .This mixture was drawn up a collection pipe, remixed, sampled, and analyzed by a nondispersive infrared analyzer.Valid time per visit and visit duration were defined as any spot measurement of CH 4 and CO 2 emissions, where the animal's head was continuously in the shroud within 20 cm of the proximity sensor for a minimum of 3 min.Total number of drops per day for each animal were also extracted.Visit data were then converted to daily emission data using SAS software program (SAS 2016) and measurements of individual animal enteric methane emissions were obtained as an average of daily CH 4 amount in gram over the test period (i.e.AVG_DAILYCH4), which was expressed as (g d −1 ).

Plasma sample collection and metabolite analyses
Near the end of GEM system tests, blood samples were collected in green topped lithium heparin vacutainers from each animal using jugular blood sampling once before morning feeding.There were no duplicates in the process of blood sample collection.Then, blood was thoroughly mixed with lithium heparin vacutainers and centrifuged to collect the plasma, which was then immediately frozen on dry ice, transported back to the lab and stored at -80˚C for subsequent metabolite analyses.For this study, blood samples of 19 heifers born in 2017 from the Angus breed, 20 heifers born in 2017 from the Charolais breed, and 20 mature cows born in 2015 from the Kinsella Composite (KC) breed were used.These animals were selected from 72, 48, and 40 cattle measured for the enteric methane emissions in 2018 to represent cattle with the lowest or highest amounts of AVG_DAILYCH4 of each breed population.

Analytical workflow
A high-performance chemical isotope labeling (CIL) method was applied to achieve relative metabolite quantification and high metabolome coverage.Fig 1 shows the overall workflow of this work.It contains the following steps: (1) sample aliquoting and generation of pooled sample, (2) 4-channel submebolome labeling, (3) LC-UV quantification and normalization, (4) mixing of 12 C-labeled individual samples and 13 C-labeled pooled sample at equal amounts, (5) LC-FTICR-MS analysis of 12 C-/ 13 C-mixtures, (6) data processing using R programs, (7) metabolite identification and statistical analysis using IsoMS Pro software (NovaMT, AB, Canada).The detailed experimental conditions are described below.

Aliquoting and generation of the pooled sample
Each individual cow plasma sample from -80˚C freezer was thawed then centrifuged at 15,000 g for 10 min.Supernatant was taken and split into 6 aliquots for four different labeling methods (15 μL each), backup (15 μL each) and preparation of pooled sample (35 μL each).For the one aliquot for preparing of pooled sample, all the aliquots from each individual sample were combined and mixed thoroughly, which was used as the reference.All of the leftovers were stored at -80˚C for future analyses.

Amine/Phenol submetabolome labeling
For the 15 μL aliquot for amine/phenol labeling, 45 μL of LC-MS grade MeOH was added to perform protein precipitation.The sample was vortexed then incubated at -20˚C for 2 hours.After that, the sample was centrifuged at 15,000 g for 10 min.45 μL of the supernatant was taken and completely dried using a Savant SC110A Speed Vac at room temperature.The sample was re-suspended in 25 μL of LC-MS grade water and ready for dansylation labeling.The labeling protocol was adapted from a previous report [8]. 25 μL of individual sample or 25 μL of the pooled sample was mixed with 12.5 μL of ACN.Then, 12.5 μL of 250 mM sodium carbonate/sodium bicarbonate buffer was added to the samples.The solution was mixed with 25 μL of freshly prepared 12 C-dansyl chloride (DnsCl) solution (18 mg/mL, in ACN) (for light labeling, individual samples) or 13 C-DnsCl solution (18 mg/mL, in ACN) (for heavy labeling, pooled sample).After incubation for 45 min at 40˚C, 5 μL of 250 mM sodium hydroxide solution was added to quench the reaction.The solution was then incubated at 40˚C for another

Carboxyl submetabolome labeling
For the 15 μL aliquot for carboxyl labeling, 45 μL of LC-MS grade ACN was added to perform protein precipitation.The sample was vortexed then incubated at -20˚C for 2 hours.After that, the sample was centrifuged at 15,000 g for 10 min.40 μL of the supernatant was taken and ready for DmPA labeling.The labeling protocol was adapted from a previous report [14].40 μL of individual sample or 40 μL of the pooled sample was mixed with 10 μL of 0.5 M triethanolamine.Then, 25 μL of freshly prepared 12 C-DmPA solution (10 mg/mL, in ACN) (for light labeling, individual samples) or 13 C-DmPA solution (10 mg/mL, in ACN) (for heavy labeling, pooled sample).After incubation for 60 min at 80˚C, 20 μL of 0.2 M Tri-Gly was added and incubated at 80˚C for another 30 min to quench the excessive labeling reagent.

Hydroxyl submetabolome labeling
For the 15 μL aliquot for hydroxyl labeling, the labeling protocol was adapted from a previous report [15].3 μl of saturated NaCl and 1.5 μL of a 6 M HCl was added to the individual sample or the pooled sample.After vortexing, 45 μL of ethyl acetate was added to extract metabolite for twice.Then, the extracted solutions were combined and completely dried.After that, the sample was re-suspended in 12.5 μL of LC-MS grade ACN.12.5 μl of freshly prepared 18 mg/ mL 12 C-(for light labeling, individual samples) or 13 C-(for heavy labeling, pooled sample) DnsCl solution (in ACN) and 12.5 μL of 24.5 mg/mL DMAP solution (in ACN) was added.The sample was vortexed, followed by spinning down, and then incubated at 60˚C for 60 minutes.After that, 2.5 μL of a 250 mM NaOH solution was added to quench the excessive labeling reagent.The solution was incubated at 60˚C for another 10 minutes.Finally, 12.5 μL of 425 mM formic acid solution in ACN/water (50:50, v/v) was added to adjust pH.

Carbonyl submetabolome labeling
For the 15 μL aliquot for carbonyl labeling, the labeling protocol was adapted from a previous report [16].45 μL of LC-MS grade MeOH was added to perform protein precipitation.The sample was vortexed, spun down, then incubated at -20˚C for 2 hours.After that, the sample was centrifuged at 15,000 g for 10 min.45 μL of the supernatant was taken and completely dried.The sample was re-suspended in 15 μL of LC-MS grade water.15 μL of 144 mM HCl solution (in MeOH) was added to the individual sample or the pooled sample.After vortexing and spinning down, 15 μL of freshly prepared 20 mM 12 C-(for the individual samples and the pooled sample) or 13 C-(for the pooled sample) DnsHz (in MeOH) was added.The sample was then vortexed, followed by spinning down.The mixture was incubated at 40˚C for 60 minutes.After that, the sample was removed from the incubator and placed in the -80˚C freezer for 10 minutes to stop the reaction.The solution was then completely dried down.Finally, the labeled metabolites were re-dissolved in 80 μL of ACN/water (50:50, v/v).

Sample mixing
According to the quantification results, 12 C-labeled individual samples and the 13 C-labeled pooled sample were mixed in equal amounts, respectively, for all the 4 channels. 12C-and 13 Clabeled pooled sample were mixed in equal amounts to serve as a quality control (QC) sample.Besides, for carboxyl, carbonyl and hydroxyl labeling, 12 C-and 13 C-labeled blanks (water) were mixed in equal volume to serve as the method blanks.The mixtures were ready for LC-MS analysis.

Data processing
After LC-FTICR-MS analysis, the entire list of centroid peaks with information of retention times, m/z values and peak intensities was exported to.csv files using Bruker Data Analysis software (Version 4.0).IsoMS [23] R-program was used to pick peak pairs, filter false-positive pairs (e.g., dimers and common adducts) and calculate peak-pair relative intensity ratios.After the alignment of peak pairs from multiple samples using the Alignment R-program, the Zerofill R-program was applied to recover the high-confidence peak pair ratios lost during the previous data processing steps.For carboxyl, carbonyl and hydroxyl labeling, Blank Subtraction R-program was used to reduce background peak pairs.

Metabolite identification
Three-tiers identification approach was used to perform metabolite identification using IsoMS Pro software.In tier 1, based on accurate mass and retention time, peak pairs were searched against a labeled metabolite library-CIL Library, which contains more than 1,300 experimental entries, including metabolites and dipeptides.In tier 2, based on accurate mass and predicted retention time, the remaining peak pairs were searched against a linked identity library (LI Library), which includes over 7,000 pathway-related metabolites and provides high-confidence putative identification results.In tier 3, based on accurate mass, the remaining peak pairs were searched against the MyCompoundID (MCID) library, which is composed of 8,021 known human endogenous metabolites (zero-reaction library) and their predicted metabolic products from one metabolic reaction (375,809 compounds) (one-reaction library) and two metabolic reactions (10,583,901 compounds) (two-reaction library).

Methane emission and metabolite data consolidation
AVG_DAILYCH4 and each metabolite concentration were examined for outliers and normality of distribution.Values that are greater or smaller than 3 times standard deviation of the mean of AVG_DAILYCH4 or metabolite concentration were excluded.Descriptive statistics of original phenotypic values of AVG_DAILYCH4 and all metabolites were provided in S1 and S2 Tables, respectively.Since animal ages at the start of measuring CH 4 and thus at blood sample collection were different, the values of AVG_DAILYCH4 and all metabolite concentrations were adjusted through a regression on the animal age factor using glm function in R programming.The regression model is shown below.
Where y is the vector of the original phenotype values, β 0 is the intercept, β 1 is the slope of the line, X 1 is the animal age in day, ε is the residual term.The residuals of the regression were used as adjusted phenotypic values for further analyses.The adjusted CH 4 (AVG_DAILYCH4) and adjusted metabolite concentrations conform to a normal distribution for each breed population as shown in normal quantile-quantile (Q-Q) plots of S1 Fig ( adjusted AVG_DAILYCH4 and adjusted metabolite A-5 as an example).

Association analyses
T-test within breed population.The animals were assigned into high or low methane emission groups within each breed based on their adjusted AVG_DAILYCH4 values.The high and low methane emissions groups had significantly different AVG_DAILYCH4 values in each of the Angus (p-value = 4.51E-05), Charolais (p-value = 2.81E-04), and KC (p-value = 3.67E-09) breed populations as shown in Table 1.
A t-test of adjusted metabolite concentrations between the high and low methane emission groups within each breed was conducted to detect metabolites that are significantly associated with the amount of methane emission.In order to consider multiple testing, permutation tests were used to determine the significance threshold value (p-value) for the t-test [24].In the permutation tests, adjusted metabolite concentrations were randomly shuffled for each metabolite and then assigned to animals in the high or low methane groups, and a t-test was conducted and the corresponding p-value was recorded.The process was repeated 1000 times and 1000 p-values were recorded.Then, all the recorded 1000 p-values were ordered from smallest to largest and the 50th p-value, which corresponds to the 5% type 1 error of the order, was used as the significance threshold value of the t-test for each corresponding metabolite.Construction of volcano plots within breed population.In order to integrate both the biological and statistical significance concerning the metabolites, volcano plots were implemented within each breed population.In statistics, a volcano plot is a type of scatter-plot that is used to quickly identify changes in large data sets composed of replicate data [25,26].In our study, volcano plots combined a measure of statistical significance (p-value) from the t-test above (between the high and low methane emission groups within each breed) with the magnitude of the change, enabling quick visual identification of those data-points (metabolites, etc.) that display large magnitude changes that are also statistically significant.The fold change (FC) of volcano plots were calculated as Mean(High group) / Mean(Low group), and the metabolites selection criteria of FC was 1 <|Log 2 FC| 6.644 (2 6.644 �100).Different from general volcano plots, we used the significance threshold value obtained from permutation tests to determine the significance of t-test in our volcano plots.The volcano plots was conducted using Metaboanalyst 5.0 and R program.
Regression analyses within breed population.A regression analysis of adjusted AVG_-DAILYCH4 to each adjusted metabolite concentration was also performed to detect metabolites that are significantly associated with the amount of methane emission.Similarly, permutation tests were used to determine the significance threshold value (p-value) at the 5% type error for the regression analysis [24].

Pathway analyses
Metabolites that were detected by both the t-test and the regression analysis as having significant associations with the enteric methane emissions in each breed population were listed for subsequent pathway analysis.Pathway analysis was conducted using Metaboanalyst 5.0 for those common detected metabolites that have a high-confidence of quantification (tier 1 and tier 2).Those commonly detected enteric methane emission associated metabolites were matched to the Bos taurus (cow) (KEGG) pathway library to identify enriched molecular processes.

Metabolite detection and identification results of four-channel labeling LC-MS
The four-channel labeling technique targets amine/phenol -, carboxyl-, hydroxyl-, and carbonyl-submetabolomes.The data was firstly processed for each submetabolome separately, and then combined together.A total of 4235 peak pairs or metabolites were detected for the 59 cow plasma samples.Three-tier identification approach was used for identifying these metabolites.In tier 1, a total of 110 peak pairs were positively identified against CIL library.In tier 2, 260 peak pairs were putatively identified with high confidence against LI library.In tier 3, 676, 1709 and 941 peak pairs were matched in the zero-, one-and two-reaction libraries, respectively, against MCID library.Thus, in total, out of 4235 detected metabolites, 3696 metabolites (87.3%) were either definitely or putatively identified.Among them, 370 metabolites were identified as high-confidence results (tier 1 and tier 2), which were more important as they can be used for other analysis (e.g., pathway analysis).The detailed metabolite detection and identification results were presented in S3 Table.

Metabolites associated with enteric methane emissions detected by the ttest and volcano plots in the three populations
Table 2 showed the number of metabolites that were significantly associated with CH 4 detected by the t-test and a complete list of metabolites names was presented in S4 Table .As summarized in Table 2, there were 1406 significant metabolites, of which 1105 were unique, from the three populations using the t-test method, with 458, 476 and 472 metabolites for Angus, Charolais and KC, respectively.And the p-value ranges were from 3.46E-05 to 9.75E-02, 7.76E-05 to 1.33E-01 and 6.12E-09 to 6.83E-02 for Angus, Charolais and KC, respectively.The volcano plots of differential metabolite identification between the high and low enteric methane emission (AVG_DAILYCH4) groups in the Angus, Charolais and KC populations were shown in Figs 2-4.As is shown in Table 2, the Log 2 FC range used for volcano plots were from -6.5918 to -1.0448 and 1.1173 to 6.6009, -6.3489 to -1.0396 and 1.0685 to 6.5368, -6.0214 to -1.0871 and 1.0063 to 6.5121 for Angus, Charolais and KC, respectively.Additional, S5 Table presented a list of the Log 2 FC, pvalues and directions for enteric methane emissions of each metabolite for the volcano plots in the Angus, Charolais, and Kinsella Composite (KC) populations.For Angus, Charolais and KC, 157, 266 and 262 differential metabolites were positively correlated with the CH 4 amount, and 265, 195 and 157 differential metabolites were negatively correlated with the CH 4 amount.Also, the number of metabolites that are significantly associated with the enteric methane emissions (AVG_DAILYCH4) detected by volcano plots were shown in Table 2, it is noticed that the volcano plots verified 422 of 458 (92.1%), 461 of 476 (96.8%) and 419 of 472 (88.8%) of the metabolites detected by the t-test in the three breed populations, which means the t-test detected additional 36, 14, and 53 metabolites in the three populations, respectively.All the metabolites identified in volcano plots were included in t-test results.As the vast majority of the metabolites detected by the t-test in the three breed populations (88.8% to 96.8%) were verified through the volcano plots and the t-test identified additional metabolites associated with CH 4 , we will focus on all the potential CH 4 associated metabolites detected by the t-test for subsequent analysis and discussion.
The number of common metabolites between different breed populations that were significantly associated with CH 4 detected by the t-test was shown in Table 3, in which 115 (14.04%) significant metabolites were common between Angus and Charolais, 149 (19.08%) significant metabolites were common between Angus and KC, and 108 (12.86%) significant metabolites were common between Charolais and KC.A Venn diagram of significant metabolites from the results of the t-test in different breed populations was shown in Fig 5, of which, 71 metabolites (6.39%) were common among the three breed populations.As is shown in S4 Table, all the 71 common metabolites had the same direction of correlations for CH 4 amount for the Angus and KC populations, of which 66 of the 71 metabolites were positively correlated with the CH 4 amount, and 5 of the 71 metabolites had the negative correlations with the CH 4 amount.Of the 71 common significant metabolites in Charolais, 66 had the negative correlations with the CH 4 amount, and 5 metabolites had the positive correlations with the CH 4 amount.However, the direction of correlations between the significant metabolites and CH 4 amount in Charolais was opposite than that in Angus and KC.

Metabolites associated with enteric methane emissions detected by regression analysis in the three populations
The number of metabolites that were significantly associated with CH 4 detected by regression was shown in Table 4 and a complete list of metabolites names was also presented in S4 Table .In total, 1651 significant metabolites, of which 1305 were unique, were discovered using the regression analysis in the three populations with 708, 478, and 465 significant metabolites in  Angus, Charolais and KC, respectively.The p value ranges in the three populations were from 1.68E-10 to 6.26E-02, 4.86E-06 to 6.49E-02 and 1.00E-08 to 5.79E-02, respectively.
As is shown in Table 3, 131 (12.42%) significant metabolites were common between Angus and Charolais, 177 (17.77%) significant metabolites were common between Angus and KC, and 94 (11.07%) significant metabolites were common between Charolais and KC.The Venn diagram (Fig 6 ) showed that 56 (4.29%) of them were common metabolites across the three breed populations.As presented in S4 Table, for Angus and KC, 49 of the common 56 metabolites were positively correlated with the CH 4 amount and 7 of the common 56 metabolites were negatively correlated with the CH 4 amount, and they had the same direction of correlations.For Charolais, 50 of the 56 common metabolites were negatively correlated with the CH 4 amount and 6 metabolites were positively correlated with the CH 4 amount.Of these 56  common metabolites in Charolais, 55 metabolites had opposite correlation directions that that in Angus and KC.The t-test method detected 265, 324, and 286 significant metabolites that were specific to the Angus, Charolais, and KC population, which represents approximately 57.9%, 68.1%, 60.6% of the significant metabolites detected in the breed population, respectively.The regression analyses discovered 456, 309, and 250 significant metabolites that were specific to the Angus, Charolais, and KC population, which represents approximately 64.4%, 64.6%, 53.8% of the significant metabolites detected in the breed population, respectively.

Comparison between the results of t-testand regression analysis
We used both the t-test and regression to detect metabolites that are significantly associated with CH 4 .As is shown in Table 5 and Fig 7, there were 394 (51%), 279 (41.3%), and 375 (66.7%) common significant metabolites between the t-test and regression analysis in the Angus, Charolais and KC populations, respectively.It was found that 64 (8.3%), 197 (29.2%), and 97 (17.3%) metabolites were detected by t-test but not by regression in the Angus, Charolais and KC populations, respectively.A total of 314 (40.7%), 199 (29.5%), and 90 (16%) metabolites were detected by regression but not by t-test in Angus, Charolais and KC, respectively.We used a high performance four-channel CIL LC-MS method in this study and were able to quantify 4235 peak pairs or metabolites using a 3-tier identification approach.We further showed the significant metabolites of high-confidence metabolite quantification (tier 1 and tier 2) detected by both the regression analysis and t-test in the three cattle populations in Table 6.As shown in Table 6, 15, 12 and 21 high-confidence detected metabolites (tier 1 and tier 2) for Angus, Charolais, and KC populations, respectively, were significantly associated with methane emissions.

Pathway analysis
A list of enriched molecular process with p-values and impact factors from the pathway analyses was presented in S6 Table for each breed population.As an illustration, plots of -log (pvalue) against the pathway impact of top four enriched molecular process were shown in S2A, S3A, and S4A Figs, respectively, for Angus, Charolais, and KC with the pathway of the highest impact and significance located on the top right corner.On the pathway schematic (S2B, S3B and S4B Figs), light blue means the metabolites cannot be identified, but it was used as a background for enrichment analysis.Those metabolites with other colors (varying from yellow to red) were the positively identified metabolites with different levels of significance.The red-colored metabolite had a more significant change between the low and high methane emission groups than a yellow-colored metabolite.The relative metabolite concentrations of those   The number of metabolites identified by t-test but not identified by regression analysis.b The number of metabolites identified by regression analysis but not identified by t-test. https://doi.org/10.1371/journal.pone.0299268.t005

Discussion
The high performance four-channel CIL LC-MS method was able to detect and quantify up to 4235 metabolites in cattle plasma samples.This high level coverage of metabolites provides a great resource to investigate biological basis of complex traits in cattle.Indeed the analyses identified more than one thousands of metabolites in total that showed significant associations with the enteric methane emissions in the three beef breed populations.However, both the results of the t-test and regression analysis revealed that methane emission associated metabolites in beef cattle were largely breed-specific.The difference of metabolites across samples was further tested by a sparse partial least squares-discriminant analysis (sPLS-DA) [27,28]that investigated distributions of samples with low CH 4 emission measurements in regarding to their metabolite concentrations in the three beef cattle population.The sPLS-DA results (S5 Fig) showed a clear separation of the three breed populations, indicating that there were significant metabolomic differences across the three breed populations.
The results of breed-specific metabolite and methane emission associations were consistent with previous reports in RNAseq studies.Robert Mukiibi et.al investigated the hepatic transcriptomic profiles and their associations with ADG, DMI, RFI and MWT in the same Angus, Charolais, and Kinsella Composite (KC) populations through global RNAseq analyses and found that the identified differentially expressed (DE) genes were largely breed-specific [29,30].Similar results were found for rumen microbiota, and some studies have also indicated that the rumen microbiota could be influenced by host breed/species in ruminants [31,32].In addition, Zhang et.al assessed the effects of breed and feed efficiency on rumen microbiota and demonstrated that cattle breed could affect rumen microbiota at both the abundance and activity level [33].
The Angus and KC populations have more common metabolites that are significantly associated with CH 4 compared with the Angus and Charolais or Charolais and KC populations although animals in the Angus and Charolais populations were all heifers and had the same diet/management levels whereas animals in the KC population were mature cows and were fed a diffident diet.It is noted that the Angus and Charolais are pure breeds while the KC population is a composite herd.Based on our breed composition analyses using DNA genotypes, 20 KC cows have 33.41% and 0.57% Angus and Charolais, respectively, plus other breeds.The greater proportion of Angus in the KC population indicates that the KC cattle are more genetically related to Angus than to Charolais, which supports the results that more significant metabolites were shared between the KC and Angus populations.Angus, a British breed, is characterized by its moderate frame and earlier fattening, which allows it to accumulate fat at an early stage, whereas Charolais is a continental European breed and is characterized by a larger frame and later maturity to fattening [34].It is reported that in comparison with the Charolais breed, Angus has greater fat depth, greater marbling score and more daily DMI but less ADG at the similar ages [18,35], presumably due to early maturity in Angus allowing the cattle to deposit more fat at a younger age [36] and in growing cattle, more energy is needed to deposit fat than protein because protein synthesis is energetically more efficient than fat synthesis [37,38].The distinctive biological processes between Angus and Charolais may explain why the directions of associations between the significant metabolites and methane emission detected in this study were largely different than that in the Angus and KC populations.
The discrepancies of the association analyses results between the t-test and the regression are likely due to the relatively small sample size used.Published reviews on the fundamental factors determining an appropriate sample size reported that sample size determination is contingent upon the significance threshold and the type of statistical test [39,40].In order to find the reason for the different results between the t-test and regression analysis in our study, the scatter plots for some uncommon significant metabolites concentration vs. CH 4 amount were presented in S6 Fig.When the relationship between the adjusted AVG_DAILYCH4 (g/per day) and metabolite concentrations is sufficiently nonlinear, which may be also due to a small sample size, the regression analysis will likely not detect the metabolite at the significance level.The small sample size may also reduce the power of detecting significant metabolites by the ttest that contrasts the means of metabolite concentrations between the low and high methane groups, but the regression analyses may detect the significant association as it regresses metabolite concentrations on methane emission amounts on all animals.
However, the small sample size should have a smaller impact in detecting metabolite and methane emissions when its association is stronger.In this study, we further examined the specific directions of metabolites and CH 4 correlations as in the S4 Table .The results suggested that the directions of all the 394, 279, and 375 common metabolite-CH 4 associations detected by the t-test are the same as that detected by the regression in Angus, Charolais and KC, respectively.Therefore, in consideration of limitations of the smaller sample size used, we further classified metabolites that were detected by both the t-test and by regression as candidate or common metabolites associated with methane emission and metabolites that were detected by one of the t-test and regression methods as suggestive.We found 394, 279 and 375 candidate or common metabolites for methane emission in Angus, Charolais and KC, and 378, 396 and 187 suggestive metabolites associated with methane emission in Angus, Charolais, and KC, respectively.
Metabolites 3-Hydroxybutyric acid (C-1373), Hydroxyisobutyric acid (C-1580) and L-Formylkynurenine/N'-Formylkynurenine (A-357) were detected in both the Angus and KC populations in Table 6.And as shown in S4 Table, 3-Hydroxybutyric acid and Hydroxyisobutyric acid had positive correlations with the CH 4 amount, and L-Formylkynurenine/N'-Formylkynurenine had a negative correlation with CH 4 amount.Minji Kim et.al [41] investigated the metabolic characteristics of Japanese Black cattle using enteric methane emissions and found that in cattle with high methane emissions the concentration of blood β-hydroxybutyric acid and insulin levels were high, whereas blood amino acid levels were low, which is similar to our results.Luz Ya ´ñez et.al found that it was feasible to produce (R)-3-hydroxybutyric acid from methane in vivo depolymerization of polyhydroxybutyrate in Methylocystis parvus OBBP [42].Mai et.al also demonstrated that 2-hydroxyisobutyric acid could be produced from methane [43].
According to the results of our study, morpholine was found to be significantly associated with the enteric methane emissions in Angus.Shaukat et.al found morpholine was a developing solvent for CO 2 capture, which had better reactivity with green gas such as CO 2 [44].In our study, phenylalanine was identified in KC and had a positive correlation with CH 4 amount.Some previous researches have shown that the phenylalanine is a promoter in methane hydrate formation [45,46].
We also listed 44 common metabolites of Tier 3 in the three beef cattle populations detected by both the t-test and regression analysis in S7 Table.These statistically significant candidate and suggestive metabolites quantified using the 3-tier approach will also help identify metabolites as biomarkers to assist with selection of cattle with reduced methane emissions.For instance, C-5074 was one of the candidate metabolites detected by both the t-test and regression analyses in all three breed populations.The histogram of relative concentration of C-5074 in high and low adjusted CH 4 group and the scatter plots for C-5074 concentration vs. CH 4 amount in three populations was shown in Fig 8 .Both of results of the t-test and regression analysis clearly indicated that the significant metabolite (C-5074) was positively associated with the amount of CH 4 in both Angus and KC but negatively associated with the amount of CH 4 in the Charolais population.Indeed, follow up studies are required to further characterize the metabolite and verify its association with methane emission through validation studies.
The metabolites detected in this study as associated with the methane emission will also help improve our understanding on the molecular basis of metabolite and methane emissions.With the help of pathway analysis, we can study and compare the importance of some pathways and metabolites, and this would facilitate our understanding of why different cows have different methane production and emission amounts.
For Angus breed (S2A Fig), the most affected pathways include "Arginine and proline metabolism", "Arginine biosynthesis", "Butanoate metabolism" and "Glutathione metabolism".We took "Arginine and proline metabolism" as an example for more discussion (S2B Fig) .L-ornithine is a non-essential, non-protein amino acid [47], and was significantly downregulated in high methane emission group.For Charolais breed (S3A Fig), "beta-Alanine metabolism", "Pyruvate metabolism", "Glycolysis / Gluconeogenesis" and "Citrate cycle (TCA cycle)" were the four most affected pathways, and as for the "beta-Alanine metabolism" (S3B Fig), we found a significant increase of 3-Oxopropanoate in high methane emission group.
For KC breed (S4A Fig), many pathways were affected and of those, "Phenylalanine, tyrosine and tryptophan biosynthesis", "Phenylalanine metabolism", "Arginine biosynthesis" and "Arginine and proline metabolism" were the four top affected ones.Here, we provide a discussion about the relation between methane emission and "Phenylalanine, tyrosine and tryptophan biosynthesis" as shown in S4B Fig. Phenylalanine is an aromatic essential amino acid.An increased concentration of phenylalanine in beef cattle with high methane emission amount was observed.A recent study also found the enrichment of phenylalanine in cows with higher milk production and less methanogen and methanogenesis functions, which may reveal the important role of the rumen microbiome [48].Arginine and proline metabolism and arginine biosynthesis were two common pathways in all four top pathway detected by the both Angus and KC populations, which also supports the results that the Angus and KC populations had more common metabolites that were significantly associated with CH 4 .It was found that the alanine pathway was the center of interactions of many other pathways [49], which implied the activation of this pathway might be related to higher methane emissions.It is reported that Arginine and proline metabolism (detected in Angus and KC population) and Phenylalanine metabolism (detected in KC population) were significantly affected pathway on enteric methane emissions from dairy cows [50].
Furthermore, the above enriched pathways or molecular processes detected through analyzing metabolite samples of the Angus, Charolais, and KC populations with methane emissions will also help prioritize candidate genes to detect DNA polymorphisms that are associated with methane emissions in beef cattle.

Conclusions
In this work, a high performance four-channel CIL LC-MS method was applied to profile amine/phenol -, carboxyl -, hydroxyl-and carbonyl submetabolomes of cow plasma samples.High metabolome coverage was achieved with the detection of 4235 metabolites in three tiers., T-tests, volcano plots, and regression analyses within breed population were conducted to detect metabolites that were significantly associated with the amount of methane emission.Our results showed that metabolites that were associated with methane emissions in beef cattle were largely breed-specific.The pathway analysis detected enriched metabolic processes related to methane production in beef cattle.The detected metabolites from the three beef cattle populations and enriched metabolic processes provide valuable resources to further characterize the metabolites and verify their associations with CH 4 as biomarkers to assist in selection of cattle with reduced methane emissions.

Fig 8 .
Fig 8.The histogram of relative concentrations of an identified metabolite (C-5074) in the high and low enteric methane emission (AVG_DAILYCH4) groups and its scatter plots of concentration vs. AVG_DAILYCH4 amount in the Angus, Charolais, and Kinsella Composite (KC) populations.https://doi.org/10.1371/journal.pone.0299268.g008

Table 1 . Significance of differences of AVG_DAILYCH4 (g/per day) between the high and low methane emission groups in the Angus, Charolais, and Kinsella Composite (KC) populations. Threshold Number of records Numb er of records in low group Number of records in high group Mean_low group±se Mean_high group±se p-value
https://doi.org/10.1371/journal.pone.0299268.t001

Table 3 . Number of common metabolites between or among the Angus, Charolais, and Kinsella Composite (KC) populations that are significantly associated with the enteric methane emissions (AVG_DAILYCH4) detected by both the t-test and regression analysis.
https://doi.org/10.1371/journal.pone.0299268.t003